
#include "../common/matrix_op.h"
#include "../common/error.h"
#include "../common/mem.h"

#include <math.h>
#include <float.h>
#include <stdlib.h>
#include <memory.h>

Vector* allocate_Vector(int n)
{
	Vector* vector;
	elem* data;

	vector = (Vector*) mem_alloc_and_zero(1, sizeof(Vector));
	data = (elem*) mem_alloc_and_zero(n, sizeof(elem));

	vector->size = n;
	vector->data = data;

	return vector;
}/* allocate_Vector */


void free_vector(Vector* vector)
{
	if (vector)
	{
		mem_free(vector->data);
		mem_free(vector);
	}/* if(vector) */

}/* free_vector */


elem scalar_mlpl(Vector* v1,Vector* v2)
{
	int i;
	elem result;
	int n;

	if (v1->size != v2->size)
		fatal_error("Trying to multiply vectors of different size!");

	n = v1->size;
	result = 0;

	for(i = 0; i < n; i++)
		result += v1->data[i] * v2->data[i];

	return result;
}/* scalar_mlpl */


void diagonal_matrix_vector_mlpl(Vector* diag,Vector* x,Vector* result)
{
	int n,i;
	n = x->size;

	for(i = 0; i < n; i++)
		result->data[i] = diag->data[i] * x->data[i];
}/* diagonal_matrix_vector_mlpl */


/* add or subtract vector a to vector b into result. NOTE: ANY VALUES IN result WILL BE OVERWRITTEN. */
/* assumes all vectors are of the same dimension and valid (not null). */
/* op: use ADD_VECTORS to add, SUBTRACT_VECTORS to subtract */
void add_vectors(Vector* a, Vector* b, Vector* result, int op)
{
	int n,i;

	n = a->size;

	for (i = 0; i < n; i++)
		result->data[i] = a->data[i] + op * b->data[i];
}/* add vectors */


elem vector_norm(Vector* vector)
{
	int n,i;
	elem norm = 0;

	n = vector->size;

	for (i = 0; i < n; i++)
	{
		norm += (vector->data[i])*(vector->data[i]);
	}

	norm = sqrt(norm);

	return norm;
}

void normalize_vector(Vector* vector)
{
	int n;
	elem norm;

	n = vector->size;

	norm = vector_norm(vector);
	if (norm != 0)
	{
		norm = 1 / norm;
		scalar_vector_mlpl(norm,vector,vector);
	}
}/* normalize_vector */

void scalar_vector_mlpl(elem scalar, Vector* x, Vector* result)
{
	int n,i;

	n = x->size;

	for (i = 0; i < n; i++)
		result->data[i] = x->data[i] * scalar;
}

void vector_copy(Vector* dstVector, Vector* srcVector)
{
	int n /*, i*/;

	n = dstVector->size;

	memcpy(dstVector->data, srcVector->data, n * sizeof(elem));

	/* for (i = 0; i < n; i++) */
	/* 	dstVector->data[i] = srcVector->data[i]; */
}

elem matrix_norm1_by_rows(elem* matrix, int size)
{
	int i, j;
	elem max, currentSum;

	max = -DBL_MAX;

	for (i = 0; i < size; i++)
	{
		currentSum = 0;

		/*  Sum the row */
		for (j = 0; j < size; j++)
			currentSum += fabs(matrix[i * size + j]);

		max = ((max > currentSum) ? max : currentSum);
	}

	return max;
}

elem matrix_norm1_by_columns(elem* matrix, int size)
{
	int i, j;
	elem max, currentSum;

	max = -DBL_MAX;

	for (j = 0; j < size; j++)
	{
		currentSum = 0;

		/*  Sum the column */
		for (i = 0; i < size; i++)
			currentSum += fabs(matrix[i * size + j]);

		max = ((max > currentSum) ? max : currentSum);
	}

	return max;
}

/* Mv */
void matrix_vector_mlpl(const elem* matrix, const Vector* vector, Vector* result)
{
	int i,j,n;
	elem temp;

	n = vector->size;

	for (i=0;i<n;i++)
	{
		temp = 0;
		for (j=0;j<n;j++)
		{
			temp += matrix[i*n+j] * vector->data[j];
		}
		result->data[i] = temp;
	}
}

int max_element(const elem* vector, const int vector_size)
{
	int i;
	elem max_value;
	int max_index;

	if (vector_size <= 0)
		return -1;

	max_index = 0;
	max_value = vector[0];

	for (i = 1; i < vector_size; i++)
	{
		if (vector[i] >= max_value)
		{
			max_index = i;
			max_value = vector[i];
		}
	}

	return max_index;
}
